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ABSTRACT 

We analyse the coarse-grained phase-space structure of the six Galaxy-scale dark 
matter haloes of the Aquarius Project using a state-of-the-art 6D substructure finder. 
Within r5o, we find that about 35% of the mass is in identifiable substructures, pre- 
dominantly tidal streams, but including about 14% in self-bound subhaloes. The slope 
of the differential substructure mass function is close to —2, which should be compared 
to ~ —1.9 for the population of self-bound subhaloes. Near r$o about 60% of the mass 
is in substructures, with about 30% in self-bound subhaloes. The inner 35 kpc of 
the highest resolution simulation has only 0.5% of its mass in self-bound subhaloes, 
but 3.3% in detected substructure, again primarily tidal streams. The densest tidal 
streams near the solar position have a 3-D mass density about 1% of the local mean, 
and populate the high velocity tail of the velocity distribution. 

Key words: methods: numerical, cosmology: dark matter 



1 INTRODUCTION 

The detailed phase-space distribution of cold dark mat- 
ter haloes can substantially affect prospects for dark mat- 
ter detection. Direct detection experiments are starting 
to probe significant fractions of the parameter space of 
plausible theoretical models, so a first detection of dark 
matter (DM) may be imminent. Only recently has re- 
alistic prediction of the phase-space structure near the 
Earth become possible, because very high-resolution nu- 
merical simulations are required. It is now well established 
that the outer regions of cold dark matter (CDM) haloes 
have a complicated phase- space structure wi t h many sub- 



haloes and tid al streams (Moore eta J 119991: iKlvpin et al 
19991 : iGhima et alJl200d:IStoehr et alJl2002b iDiemand et al 



2004. 120081 : ISprineel et alJl2008h . This raises the question 



of whether similar structures might affect direct detection 
probabilities. Could the Earth be sitting in a "hole", i.e. a 
locally very underdense region, in the DM distribution, as 
might occur in a fractal structure, or if most of the mass near 
the Sun were concentrated in dense, low-mass subhaloes. 
Some simulators have indeed argued that a significant frac- 
tion of the local mass could lie in solar or Earth-mass sub- 
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haloes fe.g. lDiemand et al] |2005), although more recent sim- 
ulations suggest that the local mass frac tion in bound sub- 
haloes of any mass is well below 1% ([Vogelsberger et all 
I2009U Vogelsberger fe WhitJ20ld ). Other possibilities might 
be for the Earth to lie within a subhalo, or within a dense 
tidal stream created by disruption of an earlier subhalo. Ei- 
ther of these would produce a spike in the velocity distribu- 
tion of local dark matter particles. 

Quantifying the detailed phase-space structure of a 
halo requires disentangling its various DM components: the 
smooth component whi ch is, in fact, a superposition o f many 
fundamental streams (|Vogelsberger fc Whitd [2010); com- 
pact, self-bound subhaloes; and tidal streams created by the 
disruption of such subhaloes. Efficient identification of the 
tidal streams requires a sensi tive and robust structure- f inder 
in 6D phase-space. Recently iMacieiewski et all ([20 09a) pre- 
sented an algorithm, Hierarchical Structure Finder (HSF), 
designed specifically for this purpose. HSF identifies struc- 
tures in an N-body simulation as coherent, overdense sets 
of particles in the full 6D position-velocity distribution. In 
this paper we use HSF to study the six Gala xy- scale halos 
simu lated as part of the Aquarius Project ([Springel et al.l 
[20Q3), providing a quantitative analysis of the various halo 
components and focussing, in particular, on the inner halo 
relevant for detection experiments. 
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Our paper is organised as follows: Section [2] presents a 
brief description of our Hierarchical Structure Finder. Sec- 
tion [3] then explores the numerical convergence of our results 
by analysing simulations of a single halo at a variety of res- 
olutions. Finally we use the full set of six Aquarius halos to 
study the expected scatter in substructure properties among 
Galaxy-scale halos. Section 3] focusses on substructure in the 
inner halo in order to assess possible consequences for direct 
detection experiments. The final section gives our conclu- 
sions. 



2 HIERARCHICAL STRUCTURE FINDER 

Many different algorithms for identifying (sub) structure 
have been applied to N-body simulations of the growth of 
cosmic structure. One of the first and most widel y used is the 
Friend s- of- Friends (FOF) scheme introduced by ID avis et al.l 
(|l985h : this defines "halos" as disjoint particle sets contain- 
ing every particle closer than some maximal linking length 
to at least one other member of its set. This bounds ob- 
jects approximately by an isodensity surface, but makes 
no assumption about their shape or internal structure. In 
contrast, the Spher ical Overdensity (SO) group-finder of 
ICole &; Lacevl ([ 19961 ) finds high-density peaks in the particle 
distribution, grows spheres centred on each until the mean 
enclosed density drops to a specified value (typically ~ 200 
times the cosmological mean) and then defines halos as the 
contents of those spheres whose centres do not lie within 
a more massive halo. More recent structure- finders [e.g. 



SKID (iGovernato et al.1 Il997h SUBF IND (ISpringel et al.1 
2001), VOBOZ (|Nevrincket al.ll2005l). PSB dKim fe PaTkl 



2006) ADAPTAHOP (lAubert Pichon fc Colombi 2004), 



AHF (iKnollmann fe Knebell2010l ). HSF (|Macieiewski et all 
2009a])] typically identify objects as connected self-bound 
particle sets above some density threshold. Such methods 
have two steps. Each particle is first linked to a local DM 
density maximum by following the gradient of a particle- 
based estimate of the underlying DM density field. The par- 
ticle set attached to a given maximum defines a candidate 
structure. In a second step, particles which are gravitation- 
ally unbound to the structure are discarded until a fully self- 
bound final object is obtained. The various methods differ 
in the way particles are treated when they belong to more 
than one candidate and in the way unbound particles are 
redistributed. Most methods produce a hierarchical charac- 
terisation of structure where halos contain subhaloes which 
in turn can contain their own subhaloes. 

These methods can be extended to higher dimensions, in 
particular to 6D phase-space. The main complication is then 
that the smoothed density field and its gradient must be 
estimated from the particle distribution in six dimensions. 
The Hierarchichal Struct ure Finder (HSF) presented by 
iMacieiewski et al . (2009a) is an algorithm of this type, and 
can be used, just like the above 3D algorithms, to identify 
bound subhaloes. Other six dimensional pha se-space struc- 
ture f inders have be en developed recently bv | Diemand et all 
(|2006l . 6dFOF) and lSharma fe Johnston! (|2009l . EnLink). In 
the following we describe the HSF method in somewhat 
more detail, since this is the method we will use for the 
rest of t his paper. Further technica l details and tests can be 
found in IMacieiewski et al.1 (|2009al ). 



To find candidate structures we first need to estimate 
phase-space densities at the positions of all the particles. 
Furthermore we need to calculate local phase-space den- 
sity gradients. HSF does this using a six-dimensional SPH 
smoothing kernel with a loca l adaptive metric as imple - 
mented in the EnBiD code ([Sharma fe Stemmed 12006). 
Neighbouring particles can then be used to derive the re- 
quired gradients. For the SPH kernel we use iV sp h = 64 
neighbours whereas for the gradient estimate we use iV n gb = 
20 neighbours. 

Once phase-space densities have been calculated, we 
sort the particles according to their density in descending 
order. Then we start to grow structures from high to low 
phase-space densities. While walking down in density we 
mark for each particle the two closest (according to the lo- 
cal phase-space metric) neighbours with higher phase-space 
density, if such particles exist. In this way we grow disjoint 
structures until we encounter a saddle point, which can be 
identified by observing the two marked particles and seeing 
if they belong to different structures. A saddle point occurs 
at the border of two structures. In the standard setup of 
HSF, which is used throughout this paper, the masses of the 
two structures separated by the saddle point are compared 
and the smaller one is cut, defining a complete individual 
structure. All particles below the saddle point whose higher 
density neighbours are part of the cut object are attached 
to the other, larger structure. Pursuing this procedure un- 
til all particles have been considered divides a halo into a 
unique disjoint set of substructures, of which the most mas- 
sive, which also contains the lowest phase-density particles, 
is the main substructure. 

When we wish to isolate self-bound subhaloes, we follow 
an identical procedure, except that each time we reach a 
saddle point, we remove all unbound particles iteratively 
from the smaller structure and attach them provisionally to 
the larger structure. Once we have followed this algorithm 
down to the lowest phase-density particle, we are left with 
a set of self-bound subhaloes and a few particles which are 
bound to no subhalo, not even the most massive self-bound 
subhalo which again is the one containing the lowest phase- 
density bound particles. 

Both these procedures divide a halo into a disjoint set of 
phase-space structures, each containing a single phase-space 
density peak and bounded approximately by a level surface 
of phase-space density. In the first procedure each structure 
normally contains both bound and unbound particles, and 
all halo particles are assigned to some structure. In the sec- 
ond procedure, each structure is self-bound, and some halo 
particles are not assigned to any structure. To be specific, 
in the following we will refer to all particles inside rscEl as 
the halo. We call the most massive substructure constructed 
from these particles the main halo. Note that by definition it 
cannot extend beyond rso and that its mean density within 
rso will be less than 50 times the critical density. Note also 
that the main halo will change slightly according to whether 
we do or do not apply the unbinding and reassignment pro- 
cedures. In the former case we refer to all other structures 



1 This is defined as the radius of the largest sphere centred on 
the halo density peak which encloses a mean density at least 50 
times the critical value. 
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as (self-bound) subhaloes, whereas in the latter case we re- 
fer to them as substructures. A subhalo is thus always part 
of a substructure, but a substructure does not necessarily 
contain a subhalo. 



3 STRUCTURES IN THE AQUARIUS 
SIMULATIONS 

We study the phase-space structure of Milky Way-sized 
DM haloes using the high-resolution simulations of the 
Aquarius Project (Spr ingel et al. 2008). The cosmological 
parameters for these simulations are Q m = 0.25, Ha = 
0.75,(78 = 0.9 and H = 73km s _1 Mpc _1 . For this project 
six Galaxy-mass haloes (Aq-A to Aq-F) were selected from 
a lower resolution version o f the Millennium-II Simulation 
(jBovlan-Kolchin et aU koOQ) and resimulated with progres- 
sively higher particle number and smaller softening length. 
The haloes were selected to have no close massive companion 
at z = 0. When studying differences in phase-space structure 
between these haloes, we use the second resolution level (the 
highest for which results are available for all six objects). At 
this resolution all haloes have more than 1.6 x 10 8 particles 
inside rso, corresponding to a particle mass ~ 1O 4 M . In 
addition, we use resimulations of the Aq-A halo at four dif- 
ferent resolution levels to check the numerical convergence 
of our results. In the final section of this paper we inves- 
tigate phase-space structure in the inner halo, defined as 
t < Tinner = 35 kpc. For this purpose we use three resolu- 
tion levels of the Aq-A halo with the largest one (Aq-A-1) 
having almost 1.5 x 10 9 particles inside rso and more than 
2 x 10 8 particles inside r^ner. Together with the ability of 
HSF to analyse the full 6D particle distribution, this simu- 
lation set allows the first robust and fully general quantifi- 
cation of the various phase-space components predicted by 
the ACDM model at r ~ 8 kpc where direct detection takes 
place. 



3.1 A resolution study 

We begin by analysing the mass functions of substructures 
and of self-bound subhaloes in the Aq-A halo and their de- 
pendence on resolution. To be consistent with earlier work 
we define the edge of the halo at r^o = 433 kpc and we 
count all objects within this radius, but we note that this is 
a large radius and, as a result, the counts are dominated by 
objects beyond 100 kpc, more than an order of magnitude 
further from the Galactic Centre than the Sun. In the upper 
panel of Fig. [T] we compare the differential mass functions 
of substructures (solid curves) and of self-bound subhaloes 
(dashed curves) at four different resolutions. The lower panel 
shows the corresponding cumulative mass functions. In both 
cases the mass functions agree quite well between the simu- 
lations above their respective resolution limits. For the self- 
bound subhaloes, the slope of the differential mass function 
is close to —1.9, as fo und earlier in the SUBFIND analysis of 
ISpringel et al.l (|20Q8h . At lower resolution the distribution is 
better approximated, particularly in the low mass bins, by 
a slope close to —1.8. For the Aq-A- 2 halo, we find that 
14% of the halo mass is in self-bound subhaloes, which is 
20% higher than the corresponding SUBFIND value (12% 
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Figure 1. Top panel: Differential substructure mass functions 
for different resolutions of the Aq-A halo in the region r < r$o = 
433 kpc. Solid lines indicate mass functions for the substructures 
identified when HSF is executed without unbinding and reassign- 
ment procedures (see the text), while dashed lines indicate the 
corresponding functions for the self-bound subhaloes found when 
these procedures are implemented. Different colours refer to the 
different resolution levels as indicated in the plot. The straight 
black line is a power-law fit, dN/dM oc M -1-9 , to the data for 
self-bound subhaloes. The data for substructures are instead fit by 
a power-law of slope —2 (the horizontal red line). Bottom panel: 
Corresponding cumulative mass functions inside r$o, expressed 
as the total fraction of the enclosed mass in identified substruc- 
ture. Line styles are the same as in the top panel. Black and red 
lines are based on the fits in the upper panel with a high mass 
truncation at 2% of the total mass. 
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- ISpringel etHl l2008h , reflecting the fact that HSF typi- 
cally attaches more mass to each identified object than SUB- 
FIND. In the lower resolution simulations, particles from un- 
resolved low-mass substructures are in many cases attached 
to more massive objects. This explains the shifts in the mass 
functions at high masses that are well resolved by all simu- 
lations (see the bottom panel of Fig[T]). The mass fraction in 
self-bound subhaloes changes from 11% for Aq-A-5 to 14% 
for Aq-A-2. 

The power-law behaviour of the mass function of self- 
bound subhaloes has been known for some time, but there 
has been controversy over its slope. If this slope is —2, then 
the mass fraction in subhaloes diverges logarithmically at 
low mass, and is cut off at a mass corresponding to the free- 
streaming length of the underlying DM particle (typically 
Earth mass for neutralino candidates). If the slope is -1.9, 
however, as found above, then the mass fraction in subhaloes 
has already effectively converged at the limit of the highest 
resolution Aquarius simulations and hence is ~ 15% within 
rso • Within smaller radii this fraction drops dramatically, as 
we will see below. 

HSF makes it possible to find unbound substructures 
also, and the solid lines in Fig.[T]show that within rso rough 
numerical convergence is achieved for their mass function. 
In this case, however, the slope appears close to —2 and the 
mass in substructures exceeds that in self-bound subhaloes 
by more than a factor of 2 at all masses, and by increas- 
ingly large amounts at small mass. To the resolution limit 
of Aq-A-2, 35% of the mass within rso is in substructure, 
showing the total mass detected in unbound tidal streams 
to be significantly larger than in self-bound subhaloes. With 
increasing resolution, significantly more substructures are 
found, and mass shifts from massive to smaller substruc- 
tures, as found above for self-bound subhaloes but even more 
strongly. The behaviour seen in the lower panel of Fig.[T]can- 
not be extrapolated straightforwardly to lower mass. As we 
will see below, at the resolution of Aq-A-2, most of the mass 
in the outer halo is resolved into substructures, but relatively 
little of the mass in the inner regions. At higher resolution it 
will be the transition between these two regimes which con- 
trols the total mass fraction in substructure, rather than the 
increase in resolved substructures at any particular radius. 

For this same resolution series of Aq-A simulations, 
Fig. [2] shows the mass fractions in substructures and in self- 
bound subhaloes for a set of of 10 disjoint spherical shells 
extending from 1 kpc to rso- Within 35 kpc we also show re- 
sults for the highest resolution simulation Aq-A-1. Although 
the results for a single simulation are rather noisy, they can 
be represented reasonably well by 



r\oc 
/sub 



: exp[7 + /31n(r/r 5 o) + 0.5a ln 2 (r/r 5 o)], 



(i) 



with parameters a — —0.31, f3 — 0.98 and 7 — —1.09 for the 
subhaloes. This analytic form was us ed to fit the radial dis - 
tribution of SUBFIND subhaloes in Spr ingel et al. I (|2QQ8h . 
We note that the quoted parameters were obtained from a 
fit to data for the full set of six level 2 Aquarius haloes (see 
below). The distribution of HSF subhaloes is very similar 
to that of SUBFIND subhaloes, but HSF attaches slightly 
more particles to objects near the centre. 

We find substructures down to 0.6 kpc in Aq-A-1 and 
down to 2 kpc in Aq-A-2 and Aq-A-3 simulations. For Aq- 
A-4, however, no substructures are found within ~ 9 kpc, 




Figure 2. Fraction of mass in substructures and in self-bound 
subhaloes as a function of radius estimated from a set of disjoint 
spherical shells and for various resolution levels of the Aq-A halo 
as indicated by colour. We take the mass within rso for the 4 
lower resolution haloes. For Aq-A-1 only the mass within the 
inner 35 kpc was used. Solid lines represent the substructures 
and dashed lines the self-bound subhaloes. The black line is an 
analytic fit based on data for the full set of 6 Aquarius haloes 
(see below) = exp[7 + ft ln(r/r 50 ) + 0.5a: In 2 (r/r 50 )] with 



parameters a = —0.31,(3 = 0.98 and 7 = 
the same function shifted vertically to 7 : 



-1.09. The red line is 
0.35. 



demonstrating that at least ~ 10 7 particles are needed 
within rso to begin to study streams around the Sun's posi- 
tion. The red curve shows the prediction of Eq. (pQ) for pa- 
rameters a = —0.31, P — 0.98 and 7 = 0.35 but, in contrast 
to the situation with subhaloes, it is clear that the results 
in the inner regions are not converging with increasing res- 
olution. The results for Aq-A-1 lie well above the red line 
and should clearly still be considered as a lower limit to the 
mass fraction contained in tidal streams in these regions. 



3.2 Mass distribution inside rso 

In this section we quantify the object-to-object scatter in the 
substructure mass function by analysing all six level 2 haloes 
of the Aquarius Project. Two of our haloes Aq-A and Aq-C 
did not experience major mergers below redshift 3. Haloes 
Aq-B and Aq-F on the other hand each underwent a major 
merger below redshift 1.5. More information on the merger 
history of the Aquarius haloes and how representative these 
haloes are of the population of Milky W ay-lik e haloes can 
be fou nd in iBovlan-Kolchin et al.l (|2009h and in I Wang et al.l 

d2oT3). 

The top panel of Fig. [3] shows the differential substruc- 
ture mass function for the six Aquarius haloes Aq-A-2 to Aq- 
F-2; the bottom panel presents corresponding cumulative 
mass functions. Table [T] lists values for the total mass frac- 
tion in substructure in each halo using different substructure 
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Halo 


SUBFIND 
(per cent) 


W Q T T "R T? T "NT "P> 


HSF subhaloes 
(per cent) 


A^HSF 


HSF substructures 
(per cent) 


Aq-A-2 


12.16 


45024 


14.14 


48052 


34.71 


Aq-B-2 


10.54 


42537 


14.44 


44143 


33.65 


Aq-C-2 


7.17 


35022 


7.72 


36525 


29.60 


Aq-D-2 


13.06 


47014 


14.26 


49726 


34.49 


Aq-E-2 


10.75 


42725 


13.10 


44400 


33.65 


Aq-F-2 


13.39 


52503 


15.16 


57269 


34.18 



Table 1. Total mass in substructures within rso for the 6 Aquarius haloes and different structure- finding methods. Mass fractions are 
calculated relative to the total mass within r$o. Halo: Aquarius halo label, SUBFIND: total mass fraction of particles in bound subhaloes 
found by SUBFIND; A^subfind-' number of bound subhaloes found by SUBFIND; HSF subhaloes: total mass fraction of particles in 
bound subhaloes found by HSF; A^hsf-' number of bound subhaloes found by HSF; HSF substructures: total mass fraction of particles 
in substructures found by HSF. 



ident ification methods. SUBFIND subhaloes ([Springel et al.l 
2008) and subhaloes found by HSF follow the same power- 
law with a slope close to —1.9. Halo Aq-C-2 evolves in a 
quiet merger environment and this explains its deficit in 
substructures. For the general substructures the slope of the 
power-law is close to —2, but it is difficult to measure this 
value accurately because the low-mass end is contaminated 
by substructure arising through discreteness noise, particles 
which are connected by HSF but do not represent physical 
substructures. 

In Aq-C-2 only 8% of the mass within rso belongs to 
self-bound subhaloes, while for Aq-E-2 this number is 13% 
and for the other haloes it is 14 — 16%. The subhalo mass 
of Aq-F-2 is dominated by the largest subhaloes. It is in- 
teresting to observe that although Aq-C-2 has the fewest 
self-bound subhaloes, its mass fraction in substructures is 
about 30%, close to the value for Aq-B-2 (34%) which had 
completely different and much richer merger history. All the 
other haloes also have substructure mass fractions around 
34%. The slope of the differential substructure mass func- 
tion is similar in all haloes except Aq-F-2, where the very 
recent merger apparently causes a bias towards massive sub- 
structures. 

Fig.|4]shows the fractions of mass in substructures (solid 
lines) and in self-bound subhaloes (dashed lines) as a func- 
tion of distance from halo centre. For most of the haloes the 
latter follows Eq. (pQ), indicated as a black solid line. Clearly, 
HSF identifies substructures close to the centre in all halos. 
This is possible because of the high density contrast in 6D 
phase-space compared to 3D configuration space. Eq. JT} 
with a different normalisation (i.e. a vertical shift, see the 
solid black line) also gives a rough fit to the radial depen- 
dence of the substructure mass fraction in most haloes, but 
we note that an independent fit of the same functional form 
(the black dashed line) suggests that the mass fraction in 
substructures increases relative to that in self-bound sub- 
haloes in the inner regions of the haloes. 

In the outskirts of the haloes, near rso, up to 30% of 
the mass is in the form of self-bound subhaloes and up to 
60% resides in substructures. 



4 THE INNER HALO 

For many years experimenters have been trying to detect 
DM in laboratory devices. Detector signals are very sensi- 



tive to the local DM phase-space distribution, so it is im- 
portant to study halo structure in detail at the solar po- 
sition. This requires high-resolution simulations like those 
of the Aquarius Project. A first phase-s pace analysis using 
the A quarius haloes was carried out by IVogelsberger et al.l 
(2009) with a focus on the local velocity and spatial distri- 
butions and their imprints on direct detection signals. Here 
we extend this study and focus on phase-space structures at 
r 2^ 8 kpc. We go bey ond the self-bound subhalo analysis of 
ISpringel et al.l (|2008f ) by using HSF, which efficiently iden- 
tifies gravitationally unbound structures like tidal streams. 
Such features can have a significant impact on DM experi- 
ments, and our goal here is to quantify the total amount of 
structure in the inner halo, both bound and unbound. We 
therefore concentrate on the region within nnner = 35 kpc 
and use the excellent resolution of Aq-A-1 to analyse struc- 
tures near the solar circle. 

4.1 Substructure in the inner halo 

To give a first impression of pha se-space structure in the in - 
ner halo we use the technique of llMacieiewski et al.l (|2009bh . 
We estimate the phas e-space density at the positi on of each 
particle with EnBiD ([Sharma &; Steinmet j [2006") and plot 
an r-v r phase-space portrait in which each pixel is colour- 
coded according to the maximum phase-space density of the 
particles it contains. The top panel of Fig. [5] shows the re- 
sulting phase-space plot for the inner part of the Aq-A-2 
halo, while the bottom panel shows the corresponding plot 
for Aq-A-1. Aq-A-2 has about 2.5 x 10 7 and Aq-A-1 about 
2 x 10 8 particles inside r in ner- The increased resolution re- 
sults in substantially more self-bound structures and tidal 
streams being visible in the inner regions of Aq-A-1. In the 
following we quantify these phase-space structures in some 
detail. 

The solid lines in Fig.[6]show the differential mass func- 
tions of the substructures found by HSF within nnner f° r the 
three highest resolution Aq-A haloes. The number of sub- 
structures more massive than 1O 7 M0 is quite small and only 
for Aq-A-1 is the dynamic range sufficient to determine a 
power- law slope, which is close to —2. While the differential 
mass function for substructures within rso (the top panel of 
Fig. [3} converges reasonably well with increasing resolution, 
this is not the case for the inner regions. Here, increasing res- 
olution enables tidal streams to be followed to significantly 
lower contrast, substantially increasing the mass attached to 
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Figure 3. Top panel: Differential subhalo mass function for the 
six different Aquarius haloes Aq-A-2 to Aq-F-2. Solid lines show 
substructures and dashed lines self-bound subhaloes. The red line 
shows a power law dN/dM oc M -1-9 which is a good fit to the 
differential mass function of the self-bound subhaloes. The black 
line shows a —2 power-law which better describes the differential 
mass function of substructures. Bottom panel: The corresponding 
cumulative mass functions. Line styles are the same as in the top 
panel. 



each one and so the total mass in substructures. This effect 
was already visible in Fig. [2] and is confirmed by the cumu- 
lative mass function shown in Fig. Here also the curves 
for different resolutions agree much less well than was the 
case when we focused on substructures within rso (see the 
bottom panel of Fig. [3}. 

Only 0.5% of the mass inside n nn er is m the form of 




Figure 4. Fraction of mass in substructure as a function of ra- 
dius for the six different haloes Aq-A-2 to Aq-F-2. Solid lines 
are for all substructures and dashed lines are for self-bound sub- 
haloes. Red lines show the mean of the halo sample and black lines 
show analytic fits using the function = exp[7 + f3 In (r/r 50) + 
0.5a ln 2 (r/r5o)] with parameters a = — 0.31,/3 = 0.98 and 7 = 
— 1.09 for subhaloes and the same function shifted vertically to 
7 = 0.35 for substructures. The black dashed line shows the best 
fit for substructures with parameters a = — 0.16, /3 = 1.10 and 
7 = 0.10. 



self-bound subhaloes. Although this number includes only 
subhaloes resolved in Aq-A-1, it is expected to increase by 
at most a factor of two i f one extrapolates d own to the 
free-streaming length (see IS pringel et al1l2008l ). The mass 
fraction in substructures in this same region increases from 
0.82% for Aq-A-3 to 3.3% for Aq-A-1. Thus tidal streams 
contain almost 7 times as much mass as self-bound sub- 
haloes at the resolution of Aq-A-1, and presumably would 
contain even more at higher resolution. As Fig. [J] shows, the 
substructure mass fraction varies as a function of radius. 
At the solar circle about 0.05% of the mass is in self-bound 
subhaloes and about 0.6% in tidal streams at the resolution 
of Aq-A-1. More than 99% of the mass appears smoothly 
distributed even at this extremely high resolution and when 
processed with a state-of-the-art 6D structure finder. 

To give a visual impression of the structures found by 
HSF, Fig. [8] shows some typical substructures in the inner 
halo of Aq-A-1. The top panel presents the main halo with 
all HSF substructures (bound and unbound) removed. In the 
second row we show the biggest bound subhalo and its at- 
tached tidal streams. These extend over nearly 60 kpc. The 
mass of this biggest substructure is 2.9 x 10 8 Mq of which 
1.9 x 10 8 M is in the self-bound subhalo. In the middle 
row we show an example of the same kind of substructure, 
but here the tidal tails are more pronounced and the stream 
passes within 8 kpc of halo centre. The bottom two rows 
show two typical stream-like structures which have no asso- 
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Figure 5. Phase-space r-v r plots of Aq-A-2 (top) and Aq-A-1 
(bottom). These were constructed by first calculating the phase- 
space density at each particle using EnBiD ( Shar ma &; Steinmetd 
2006). We then create a pixelised image with 1000 x 1000 bins 
and for each pixel we calculate the logarithm of the maximum 
phase-space density over all particles in that pixel. The phase- 
space density is measured in units of Mq kpc -3 km -3 s 3 . 
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Figure 7. Cumulative mass functions for the inner part of the Aq- 
A halo (r < ri nner = 35 kpc) at different resolutions. Solid lines 
represent substructures and dashed lines self-bound subhaloes. 
All curves are normalised to Mj nner , the mass within rj nner . 



ciated self-bound subhalo. Their masses are similar to those 
of the biggest subhaloes. 
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Figure 8. Plots of the main halo and some selected substructures 
in x-y and r-v r projections for Aq-A-1. We create a 500 x 500 
image and colour each pixel according to the logarithm of the 
maximum 6D phase-space density over the enclosed particles as 
estimated using EnBiD. The phase-space density is measured in 
units of Mq kpc -3 km -3 s 3 . The mass of each substructure is 
given in units of Mq at the bottom right of each row. 



4.2 Density probability distribution function 

In this section we try to understand how the different com- 
ponents contribute to the dark matter density near the Sun. 
We ask how likely is it for the Solar System to lie within a 
subhalo or tidal stream of given local density. This is accom- 
plished by computing probability density distributions for 
the space density at random points with in a thick spherical 
shell c entred at the Solar radius, as in IVogelsberger et a"D 
(2009), but separating the dark matter into the different 
phase-space components identified by HSF. The result of 
this procedure is presented in Fig. [9] 

The left panel of Fig. [9] shows the probability distribu- 
tion function of DM density for various structures in the 
inner halo. To make this plot we first estimated the density 
at the position of every particle with radius between 6 and 
12 kpc using a standa rd SPH scheme based on 32 neigh- 
bours. As described in Vogelsb erger et ahl ((2009), we then 
fitted a smooth model to these values assuming the density 
to be stratified on similar, concentric ellipsoids and to be a 
power law of radius. This defines a model density p s heii at 
the position of each particle which can be compared with the 
directly estimated local density. This step is crucial to ac- 
count for the large density gradients in the inner halo so that 
we can focus on small-scale variations due to substructure. 
We then repeat the SPH density estimates for the subsets of 
particles in this radial range corresponding to each individ- 
ual subcomponent: the main subhalo and each individual 
self-bound subhalo and substructure. In the following we 
plot all density distributions as functions of p/p s heii and we 
construct volume-weighted probability distributions by his- 
togramming the particles with individual weights m p /pV, 
where m p is the particle mass and V the total volume be- 
tween 6 and 12 kpc. For the total mass distribution and the 
main subhalo the resulting distributions give the probability 
that an observer at a random point in this radial range will 
see local density contrast p/p s he\\- For the self-bound sub- 
halo and substructure components, the distributions show 
the mean number of self-bound subhaloes or substructures 
with local density contrast p/p s he\\ at a random point. 

If we consider first the probability distribution of den- 
sity contrast for the total mass (i.e. the sum of all the compo- 
nents) we see that there is a lognormal distribution centred 
on p/psheii = 1 with an additional low amplitude, power- 
law tail to high densities . This result was already given in 
IVogelsberger et al.1 (2009). As they showed, the lognormal 
part of the distribution reflects discreteness noise in our den- 
sity estimator. We demonstrate this again here by plotting 
as a dashed black line the distribution of analogous den- 
sity estimates for points sampled from a uniform Poisson 
distribution. This line cannot be distinguished from that 
corresponding to the main subhalo in Fig. [9j demonstrat- 
ing that this component follows our ellipsoidal model very 
closely. The power law tail is due to self-bound subhaloes, 
as is evident from its close agreement with the distribution 
calculated directly from this component. This general be- 
haviour was pointe d out not only in t he an alysis of this 
same simulation by Vogel sberger et all (l2009fh but also in 
the analytic model of Kamionkowski &; Koushiappad ([2008) 
an d in the analysis of a differe nt high resolution simulation 
by iKamionkowski et. al.l (|2010h . 

If we now consider the density contrast distributions 
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Figure 9. Left panel: Volume-weighted density probability distributions for particles in the radial range 6 — 12 kpc in Aq-A-1. For each 
particle an SPH-density is calculated using 32 neighbours. The resulting density field is fitted to a smooth ellipsoidal power-law model 
to obtain p s h e ii- SPH-densities are also calculated using the particles in each individual subhalo and substructure separately; we do not 
consider individual subcomponents containing fewer than 64 particles. The results can be used to derive a density contrast p/p s he\\ at 
each particle's position, both in total mass and for the individual subcomponent to which the particle belongs. If V denotes the total 
volume between 6 and 12 kpc and m p is the particle mass, then m p /pV is the probability that a random point in this radial range 
overlaps the particle. By histogramming this quantity for all particles we obtain the probability that a random point has density contrast 
p/pshell (the black solid line). By instead histogramming this quantity for the particles in a single component using the p values calculated 
for individual subcomponents we obtain the mean number of subcomponents at a random position with local density contrast p/p s hell 
(yellow solid curve for the main subhalo, green for the set of self-bound subhaloes, red for the set of all substructures). Labels give the 
total mass in each component (in units of Mq). The black dashed line indicates the density contrast distribution produced by our density 
estimator for a Poisson realisation of a uniform density field. Right panel: Fraction of the mass in the radial range 6 to 12 kpc with 
density contrast above p/p s hell f° r all particles, for the main subhalo, for the self-bound subhalo population and for the substructure 
population. For the latter two, the density contrast is that of the individual object containing the particle. Colours are as in the left 
panel. This panel also shows cumulative mass fraction plots for the five most massive self-bound subhaloes (long-dashed lines) and for 
the five most massive substructures (dotted lines). The masses of the individual objects are given in parentheses (in units of Mq). 



of the individual components, we see that self-bound sub- 
haloes are detectable not only in the high-density tail but 
also down to contrasts as small as 10 -4 . This reflects the 
excellent resolution of the Aq-A-1 simulation and, more im- 
portantly, the fact that our 6-D structure finder can iden- 
tify subhalo material even at very low density contrast be- 
cause of its small internal velocity dispersion. HSF identifies 
general phase-space substructure (e.g. tidal streams) down 
to even lower contrasts, of order 10 -7 . The additional sub- 
structure mass which is not part of self-bound subhaloes is 
almost entirely in this low-contrast regime. It is interesting 
that the "probability" a random point lies in such a low- 
density tidal stream reaches values much larger than one, 
meaning that HSF has identified multiple structures at each 
point in 3-space. Note, however, that because low-density 
tidal streams have an effective spatial dimensionality less 
than 3, their p/p s heii values are biased low (and the mean 
stream number correspondingly high) by the spherical kernel 
of the SPH density estimator. In the radial range between 
6 and 12 kpc selected for this analysis about 0.09% of the 
mass is in the form of self-bound subhaloes and about 0.7% 
in substructure. Thus low-density tidal streams account for 
almost 90% of the substructure detected by HSF. 

The right panel of Fig. [9] shows, for the various compo- 



nents, a cumulative plot of the mass at local density con- 
trast exceeding p/p s heii, expressed as a fraction of the total 
mass between 6 and 12 kpc. Here we see explicitly that the 
substructure component contains almost ten times as much 
mass as the bound subhalo component and that this excess 
lies almost exclusively at contrasts below 0.1. This panel 
also gives similar cumulative data for the five individually 
most massive self-bound subhaloes and for the five individu- 
ally most massive substructures. Almost a third of the mass 
in self-bound subhaloes is contained in these five objects, 
almost all of it at density contrasts exceeding unity. How- 
ever, only the most massive subhalo corresponds to one of 
the five most massive substructures, accounting for most of 
its mass. The other four massive substructures are unbound 
tidal streams with no associated subhalo. The maximum 
density contrast of these unbound streams is ~ 10 -2 (again 
this is probably biased low). The five most massive sub- 
structures together account for only about 10% of the total 
substructure mass. 

4.3 Velocity distributions 

Not only the mass density, but also the velocity distribu- 
tion of DM particles in the vicinity of the Earth is relevant 
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for direct detection experiments. Vogelsb erger et all (2009) 
showed that although this velocity distribution is quite well 
approximated by a smooth trivariate Gaussian, potentially 
measurable features are imprinted on the corresponding en- 
ergy distribution by the detailed formation history of the 
Milky Way's halo. Here we concentrate on the velocities of 
the different phase-space components in the inner halo. In 
Fig. [10] we show v r -v t projections of the distribution of all 
particles in the radial range 6 to 12 kpc (top), of those in 
substructures (middle), and of those in self-bound subhaloes 
(bottom). These plots are two-dimensional histograms, with 
colour encoding the mass in the corresponding bin as indi- 
cated by the colour bar (in units of solar masses per 2 km/s 
x 2 km/s pixel). The total mass contributing to each panel 
is given in its top right-hand corner. The main halo and 
so the bulk of the particles lie primarily at velocities below 
200 km s _1 , whereas subhaloes and tidal streams are found 
almost exclusively at higher velocities. As a result, the most 
massive subhaloes are still (just) visible in the top panel 
despite the fact that they contribute less than a tenth of 
a percent of the mass. These structures contribute to the 
high-energy tail of the recoil spectrum in direct DM detec- 
tion experiments, and so may be visible in high resolution 
experiments, particularly those with directional sensitivity 
which can detect the common motion of the substructure 
particles. 



5 CONCLUSIONS 

We study the population of subhaloes and tidal streams 
in six Milky Way-like DM haloes taken from the Aquar- 
ius Project. These structures ar e identified using the Hier - 
archical Structure Finder fHSF [Maciejewski et al.l l2009ah , 
a state-of-the-art structure finder which operates in 6-D 
phase-space. 

We find that that the differential mass function of self- 
bound subhaloes can be well described by a power-law with 
slope close to —1.9. This agrees with results from an inde- 
pendent analysis u sing the 3-D structure finder SUBFIND 
(jSpringel et al.ll2008h . Typically HSF attaches slightly more 
particles to subhaloes than SUBFIND, and also finds slightly 
more subhaloes above the simulation resolution limit (see 
Table [J). This agrees w ith previous results described in 
iMacieiewski et alj (|2009al ). About 14% of the mass within 
r5o is in self-bound subhaloes, with significant scatter among 
the six haloes Aq-A to Aq-F. HSF subhalo masses are ~ 10% 
larger than those found by SUBFIND, although the increase 
can be larger near halo centre. In most haloes the total sub- 
halo mass is dominated by the largest objects. The radial 
distributions of HSF and SUBFIND subhaloes are almost 
identical, although HSF can identify subhaloes closer to halo 
centre due to their enhanced density contrast in phase-space. 

The differential mass function for substructures (i.e. 
both subhaloes and tidal streams) is also well described with 
a power- law, but in this case the slope is close to —2. This 
is independent of simulation resolution and holds approxi- 
mately for all six haloes, thus appearing robust. For most of 
the level 2 haloes around 35% of the mass within rso is as- 
signed to substructures with mass above - 3x 1O 5 M . The 
radial distribution of these objects can be appro ximated by 
equation (pQ) introduced in lSpringel et al.l (|2QQ8h but with a 
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Figure 10. Radial velocity v r - tangential velocity vt phase-space 
plots for the 6 — 12 kpc region of the Aq-A-1 halo. For each plot 
a two-dimensional histogram was calculated and colours were set 
to reflect the mass in each 2 km s _1 x 2 km s — 1 pixel as shown 
by the colour bar, labelled in units of Mq. The top, middle and 
bottom panels show the distributions for all particles, for sub- 
structures and for subhaloes respectively. The low velocity region 
is dominated by the main halo but some substructures/subhaloes 
are still visible in the the high v r and vt velocity regions. HSF 
only detects significant substructure in these regions of phase- 
space. The total component mass in Mq is indicated in the top 
right-hand corner of each panel. 
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higher normalisation than applies for self-bound subhaloes. 
In the inner halo almost 10 times as much mass is detected 
in unbound tidal streams as in self-bound subhaloes. This 
reflects the efficiency with which tidal forces destroy bound 
subhaloes in the inner halo. 

In our highest resolution halo, Aq-A-1, HSF assigns 
about 0.5% of the mass within 35 kpc to self-bound sub- 
haloes, and about 3.3% to substructures (subhaloes and 
tidal streams), with masses higher then 3 x 10 4 Mq. In this 
region the largest phase-space substructures are either self- 
bound subhaloes with massive tidal tails stretching across 
the entire inner 35 kpc region or equally massive tidal 
streams with no attached subhalo at r < 35 kpc. The largest 
individual substructures in the inner region have masses up 
to 3 x 10 8 M^. 

IVogelsberger et al.1 (2009) showed that the density field 
in the radial range from 6 to 12 kpc is very well represented 
by a simple smooth model where density is stratified on 
similar concentric ellipsoids and falls as a power law of ra- 
dius. Fluctuations around this model are small except for a 
low-amplitude power-law tail to high density contrast which 
corresponds to self-bound subhaloes. Our HSF analysis is 
consistent with these results and allows us, in addition, to 
study the contrast of individual substructures with respect 
to the smooth background, the main subhalo, in which they 
are embedded. HSF is able to identify not only the high- 
contrast cores of individual self-bound subhaloes, but also 
their outskirts where the density contrast drops to values 
as low as - 10~ 4 . The maximum contrast of unbound tidal 
streams is ~ 10 -2 . Since these streams contain 0.6% of the 
total mass between 6 and 12 kpc, of order one massive tidal 
stream is predicted to pass through every point, contributing 
a few tenths of a percent of the local DM density. In con- 
trast, only 0.09% of the mass in this region is contributed 
by self-bound subhaloes, and the chance that the Earth lies 
in the high-density contrast region of such a subhalo is be- 
low 10 -4 . Both subhaloes and tidal streams populate the 
high-energy tail of the velocity distribution preferentially, 
and would show up in direct DM detection experiments as a 
small but significant part of the signal in events with almost 
identical (vector) momenta. 
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